function PlotXGlomOutputDist_version3(voxels,diams,DatasetID,PixelWidth,ThresholdRange, DiamsRange,MRorCT, GlomStats)
%% Plot of the statistical output of the Xglom pipeline
%
%
%   - Michael Eager,   (michael.eager@monash.edu) 
%   - (c) 2012, Monash Biomedical Imaging, Monash University, Australia

%     Copyright © 2012-2013 Michael Eager <michael.eager@monash.edu> 
%
%     This file is part of Xglom.
% 
%     This is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

if nargin < 4
 PixelWidth=0.0508; %mm, default for MR Agilent 9.4T,   Synch CT resolution 0.0153 mm 
end



if nargin < 5
    ThresholdRange = 0.11:0.01:0.5;
end
if nargin < 6
    if PixelWidth > 0.04
    DiamsRange = 1:0.1:15;
    else 
    DiamsRange =1:0.1:300;
    end
end
if nargin < 7
    MRorCT=0;
end
if nargin < 8
  GlomStats = AnalyseThresholdStatistics(diams, voxels, PixelWidth, ThresholdRange, DiamsRange,MRorCT)
end

figure1 = figure;
set (figure1, 'position',[98    52   630   830])

% Grey lines
numthresh = size(voxels,1);
greycolor= 0.6:0.2/(numthresh-1):0.8
greycolor= [greycolor' greycolor' greycolor']

%% Subplot 1: Histogram of voxel cluster sizes
subplot(2,1,1);
%stairs((voxels'),'LineWidth',0.5,'Color',[0.8 0.8 0.8]);
%hold('on')
stairs(voxels(GlomStats.optimal_index,:),'-b','LineWidth',2);
x=1:(size(voxels,2)-1);
g1 =  (GlomStats.GlomVoxelOptimum.fitobj.a1).*exp(-((x-GlomStats.GlomVoxelOptimum.fitobj.b1)./(GlomStats.GlomVoxelOptimum.fitobj.c1)).^2);
    g2 =  (GlomStats.GlomVoxelOptimum.fitobj.a2).*exp(-((x-GlomStats.GlomVoxelOptimum.fitobj.b2)./(GlomStats.GlomVoxelOptimum.fitobj.c2)).^2);
    g3 =  (GlomStats.GlomVoxelOptimum.fitobj.a3).*exp(-((x-GlomStats.GlomVoxelOptimum.fitobj.b3)./(GlomStats.GlomVoxelOptimum.fitobj.c3)).^2);
%    plot(x,g1,'-g',x,g2,'-b',x,g3,'-r', x,g1+g2+g3,'-k')
hold('off')
title(['Xglom --  Glom Size & Diameter Histograms : ' DatasetID],'FontSize',18)
if PixelWidth < 0.04
  xlim([0 100]);ylim([0 3000])
else
  xlim([0 60])  
vox_mean_height=ceil(GlomStats.GlomVoxelOptimum.mean);
if vox_mean_height ~=0
maxylim = 5*voxels(GlomStats.optimal_index,vox_mean_height)
if maxylim ==0
    maxylim=5000
end
ylim([0 2000])

end
end
xlabel('(Glomeruli) Cluster Sizes (# voxels)','FontSize',14)
ylabel('Number of (Glomeruli) Clusters','FontSize',14)

% annotation('textbox',...
% [0.581357142857143 0.633333333333333 0.277571428571429 0.102380952380957],...
%     'String',{'Estimated Glomeruli count:', num2str(sum((voxels(25,:)')))},...
%     'FitBoxToText','off');
 % imagesc(nimage)
% colormap(flipud(gray));
%axis square;

% %% Subplot 2:  Histogram of equivalent diameter
% subplot(3,1,2);
% stairs(DiamsRange*PixelWidth*1000,(diams'),'LineWidth',0.5,'Color',[0.8 0.8 0.8]);hold ('on')
% stairs(DiamsRange*PixelWidth*1000,diams(GlomStats.optimal_index,:),'-b','LineWidth',2);
% x=DiamsRange*PixelWidth*1000
% g1 =  (GlomStats.GlomVoxelOptimum.fitobj.a1).*exp(-((x-GlomStats.GlomDiamOptimum.fitobj.b1)./(GlomStats.GlomDiamOptimum.fitobj.c1)).^2);
%     g2 =  (GlomStats.GlomDiamOptimum.fitobj.a2).*exp(-((x-GlomStats.GlomDiamOptimum.fitobj.b2)./(GlomStats.GlomDiamOptimum.fitobj.c2)).^2);
%     g3 =  (GlomStats.GlomDiamOptimum.fitobj.a3).*exp(-((x-GlomStats.GlomDiamOptimum.fitobj.b3)./(GlomStats.GlomDiamOptimum.fitobj.c3)).^2);
%     plot(x,g1,'-g',x,g2,'-b',x,g3,'-r', x,g1+g2+g3,'-k')
% 
% 
% if PixelWidth < 0.04
%     xlim([0 300]);ylim([0 3000])
% else
%     xlim([DiamsRange(1) DiamsRange(end)]*PixelWidth*1000)
%     %[75 300]); %ylim([0 10000])
%     GlomStats.GlomDiamOptimum.histrange(2)
%     diam_mean_height=ceil(GlomStats.GlomDiamOptimum.mean);
%     if diam_mean_height ~=0
%         maxylim = 3*voxels(GlomStats.optimal_index,diam_mean_height)
%         if maxylim ==0
%             maxylim=5000
%         end
%         
%     end
%     %ylim([0 5000])
% end
% 
% 
% hold('off')
% title (['Glomeruli  Diameter Histogram  (Voxel Resolution ' num2str(PixelWidth*1000) ' um)'],'FontSize',14)
% ylabel('Number of Clusters','FontSize',14)
% xlabel('Eqivalent Diameter (Spherical Volume Estimate) (microns) ','FontSize',14)


%% Subplot 3
subplot(2,1,2);%().*PixelWidth*1000
if length(ThresholdRange) == size(sum((diams')),1)
    plot(ThresholdRange,sum((diams')))
else
plot(GlomStats.areas_estimate) %sum((diams')))
end
%xlim([0 300]); ylim([0 10000])
title ('Optimal Thresholding Curve','FontSize',14)
ylabel('Number of Clusters','FontSize',14);
xlabel('Threshold','FontSize',14);
yl = ylim; ymax=yl(2);
pos = get(gca,'Position');




%% Text boxes
ci = confint(GlomStats.GlomVoxelOptimum.fitobj,0.95)
ind=GlomStats.GlomVoxelOptimum.varindx;
if ~ind
    ind=2;
end
vstr=  {'Voxel Estimation Method: ',...
    [' Glomeruli Count :  ' num2str(ceil(GlomStats.GlomVoxelOptimum.total), '%d')],...
    [' Histogram Range :  ' num2str(GlomStats.GlomVoxelOptimum.histrange) ],...
    [' Voxel Size Mean  (95% CI) :  ' num2str(GlomStats.GlomVoxelOptimum.mean, '%2.2f') ' (' num2str(ci(1,ind), '%2.2f') ',' num2str(ci(2,ind), '%2.2f') ')' ],...
    [' Voxel Size StDev (95% CI) :  ' num2str(GlomStats.GlomVoxelOptimum.sd, '%2.2f')  ' (' num2str(ci(1,ind+1), '%2.2f') ',' num2str(ci(2,ind+1), '%2.2f') ')' ]}
%    [' Eq. Volume Mean  (95% CI) :  ' num2str(GlomStats.GlomVoxelOptimum.mean * (PixelWidth)^3, '%.4g') ' (' num2str(ci(1,ind) * (PixelWidth)^3, '%.4g') ',' num2str(ci(2,ind) * (PixelWidth)^3, '%.4g') ')' ],...
%    [' Eq. Volume StDev (95% CI) :  ' num2str(GlomStats.GlomVoxelOptimum.sd * (PixelWidth)^3, '%.4g')  ' (' num2str(ci(1,ind+1) * (PixelWidth)^3, '%.4g') ',' num2str(ci(2,ind+1) * (PixelWidth)^3, '%.4g') ')' ]};
% Create textbox
annotation(figure1,'textbox',...
    [0.45 0.802995391705069 0.6 0.0967741935483871],...
    'String',vstr,...
    'FitBoxToText','off', 'edgecolor','none',...
    'FontSize',12);


ci = confint(GlomStats.GlomDiamOptimum.fitobj,0.95)
ind=GlomStats.GlomDiamOptimum.varindx;
if ~ind
    ind=2;
end
dstr =  {['Diameter Estimation Method: ']...
    [' Glomeruli Count :  ' num2str(ceil(GlomStats.GlomDiamOptimum.total), '%d')],...
    [' Histogram Range :  ' num2str(GlomStats.GlomDiamOptimum.histrange)],...
    [' Diameter Mean (95% CI): ' num2str(GlomStats.GlomDiamOptimum.mean, '%2.2f')  ' (' num2str(ci(1,ind), '%2.2f') ',' num2str(ci(2,ind), '%2.2f') ')' ],...
    [' Diameter StDev (95% CI): ' num2str(GlomStats.GlomDiamOptimum.sd, '%2.2f')  ' (' num2str(ci(1,ind+1), '%2.2f') ',' num2str(ci(2,ind+1), '%2.2f') ')' ]};
% Create textbox
% annotation(figure1,'textbox',...
%     [0.45 0.505760368663594 0.6 0.0944700460829496],...
%     'String',dstr,...
%     'FitBoxToText','off',...
%     'EdgeColor','none',...
%     'FontSize',12);

disp(GlomStats.GlomVoxelOptimum.fitobj)
%disp(GlomStats.GlomDiamOptimum.fitobj)



%% Peak values and cumulative histogram
% 
% optimal_thresh = ThresholdRange(GlomStats.optimal_index);
% 
% [hmax imax] = hist(tmpimg(:),300);
% chmax=cumsum(hmax(end:-1:1));
% figure;
% [AX,H1,H2]=plotyy(imax(2:end),hmax(2:end),imax(end:-1:2),chmax(1:end-1));
% set(get(AX(1),'Ylabel'),'String','Number of Gloms');
% set(get(AX(2),'Ylabel'),'String','Cumulative Histogram of Gloms');
% xlabel('Normalised Magnitude of Peak') 
% title('Glom Peak Distribution and Cumulative Histogram (Reverse)','FontSize',14) 
% 
